arXiv:l 502.01988v2 [math.ST] 21 Oct 2015 


Computational and Statistical Boundaries for Submatrix 
Localization in a Large Noisy Matrix 

T. Tony Cay* Tengyuan Liang ' and Alexander Rakhl i n - 


Department of Statistics 
The Wharton School 
University of Pennsylvania 


Abstract 

The interplay between computational efficiency and statistical accuracy in high-dimensional in¬ 
ference has drawn increasing attention in the literature. In this paper, we study computational and 
statistical boundaries for submatrix localization. Given one observation of (one or multiple non¬ 
overlapping) signal submatrix (of magnitude A and size k m x k n ) contaminated with a noise matrix 
(of size m x n ), we establish two transition thresholds for the signal to noise A/cr ratio in terms of 
m, n, k m , and k n . The first threshold, SNR C , corresponds to the computational boundary. Below 
this threshold, it is shown that no polynomial time algorithm can succeed in identifying the subma¬ 
trix, under the hidden clique hypothesis. We introduce adaptive linear time spectral algorithms that 
identify the submatrix with high probability when the signal strength is above the threshold SNR C . 
The second threshold, S N R s , captures the statistical boundary, below which no method can succeed 
with probability going to one in the minimax sense. The exhaustive search method successfully finds 
the submatrix above this threshold. The results show an interesting phenomenon that SNR C is al¬ 
ways significantly larger than SNR s , which implies an essential gap between statistical optimality 
and computational efficiency for submatrix localization. 


1 Introduction 

The “signal + noise" model 


X = M+Z, (1) 

where M is the signal of interest and Z is noise, is ubiquitous in statistics and is used in a wide range 
of applications. When M and Z are matrices, many interesting problems arise under a variety of struc¬ 
tural assumptions on M and the distribution of Z. Examples include sparse principal component anal¬ 
ysis (PCA) (Vu and Lei, 2012; Berthet and Rigollet, 2013b; Birnbaum et al., 2013; Cai et al., 2013, 2015), 
non-negative matrix factorization (Lee and Seung, 2001), non-negative PCA (Zass and Shashua, 2006; 
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Montanari and Richard, 2014). Under the conventional statistical framework, one is looking for optimal 
statistical procedures for recovering the signal or detecting its presence. 

As the dimensionality of the data becomes large, the computational concerns associated with statis¬ 
tical procedures come to the forefront. In particular, problems with a combinatorial structure or non- 
convex constraints pose a significant computational challenge because naive methods based on exhaus¬ 
tive search are typically not computationally efficient. Trade-off between computational efficiency and 
statistical accuracy in high-dimensional inference has drawn increasing attention in the literature. In 
particular, Chandrasekaran et al. (2012) and Wainwright (2014) considered a general class of linear in¬ 
verse problems, with different emphasis on convex geometry and decomposition of statistical and com¬ 
putational errors. Chandrasekaran and Jordan (2013) studied an approach for trading off computational 
demands with statistical accuracy via relaxation hierarchies. Berthet and Rigollet (2013a); Ma and Wu 
(2013); Zhang et al. (2014) focused on computational requirements for various statistical problems, such 
as detection and regression. 

In the present paper, we study the interplay between computational efficiency and statistical accu¬ 
racy in submatrix localization based on a noisy observation of a large matrix. The problem considered 
in this paper is formalized as follows. 

1.1 Problem Formulation 

Consider the matrix X of the form 


X = M+Z, where M = A-Ir l£ (2) 

and lR m e K'" with 1 on the index set R m and zero otherwise. Here, the entries Zjj of the noise matrix 
are i.i.d. zero-mean sub-Gaussian random variables with parameter a (defined formally in Equation 
(4)). Given the parameters m, n, k m , k n ,A/a, the set of all distributions described above—for all possible 
choices of R m and C n —forms the submatrix model Mini, n, k m , k n , A la). 

This model can be further extended to the case of multiple submatrices as 

r 

m=Xa s -i Rs i£ 5 (3) 

5=1 


where \R S \ = k[ m] and C s \ = Ic" 1 denote the support set of the 5-th submatrix. For simplicity, we first 
focus on the single submatrix and then extend the analysis to the model (3) in Section 2.5. 

There are two fundamental questions associated with the submatrix model (2). One is the detection 
problem: given one observation of the X matrix, decide whether it is generated from a distribution in 
the submatrix model or from the pure noise model. Precisely, the detection problem considers testing of 
the hypotheses 

Hq:M = 0 v.s. H a :Me^(_m,n,k m ,k n ,A/a). 

The other is the localization problem, where the goal is to exactly recover the signal index sets R m and 
C n (the support of the mean matrix M). It is clear that the localization problem is at least as hard (both 
computationally and statistically) as the detection problem. As we show in this paper, the localization 
problem requires larger signal to noise ratio Ala, as well as a more detailed exploitation of the submatrix 
structure. 

If the signal to noise ratio is sufficiently large, it is computationally easy to localize the submatrix. 
On the other hand, if this ratio is small, the localization problem is statistically impossible. To quantify 
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this phenomenon, we identify two distinct thresholds (SNR s and SNR C ) for Alo in terms of parameters 
m, n, k m , lc „. The first threshold, SNR s , captures the statistical boundary, below which no method (pos¬ 
sibly exponential time) can succeed with probability going to one in the minimax sense. The exhaus¬ 
tive search method successfully finds the submatrix above this threshold. The second threshold, SNR C , 
corresponds to the computational boundary, above which an adaptive (with respect to the parameters) 
linear time spectral algorithm finds the signal. Below this threshold, no polynomial time algorithm can 
succeed, under the hidden clique hypothesis, described later. 

1.2 Prior Work 

There is a growing body of work in statistical literature on submatrix problems. Shabalin et al. (2009) pro¬ 
vided a fast iterative maximization algorithm to solve the submatrix localization problem. However, as 
with many EM type algorithms, the theoretical result is very sensitive to initialization. Arias-Castro et al. 
(2011) studied the detection problem for a cluster inside a large matrix. Butucea and Ingster (2013); Bu- 
tucea et al. (2013) formulated the submatrix detection and localization problems under Gaussian noise 
and determined sharp statistical transition boundaries. For the detection problem, Ma and Wu (2013) 
provided a computational lower bound result under the assumption that hidden clique detection is com¬ 
putationally difficult. 

Balakrishnan et al. (2011); Kolar et al. (2011) focused on statistical and computational trade-offs for 
the submatrix localization problem. They provided a computationally feasible entry-wise thresholding 
algorithm, a row/column averaging algorithm, and a convex relaxation for sparse SVD to investigate the 
minimum signai to noise ratio that is required in order to localize the submatrix. Under the sparse regime 
k m ^ m 1/2 and k n ^ n 112 , the entry-wise thresholding turns out to be the “near optimal” polynomial¬ 
time algorithm (which we will show a de-noised spectral algorithm that perform slightly better in Sec¬ 
tion 2.4). However, for the dense regime when k m ^ m 1/2 and k n £3 n 112 , the algorithms provided in Kolar 
et al. ( 2011 ) are not optimal in the sense that there are other polynomial-time algorithm that can succeed 
in Ending the submatrix with smaller SNR. Concurrently with our work, Chen and Xu (2014) provided a 
convex relaxation algorithm that improves the SNR boundary of Kolar et al. (2011) in the dense regime. 
On the downside, the implementation of the method requires a full SVD on each iteration, and there¬ 
fore does not scale well with the dimensionality of the problem. Furthermore, there is no computational 
lower bound in the literature to guarantee the optimality of the SNR boundary achieved in Chen and Xu 
(2014). 

A problem similar to submatrix localization is that of clique finding. Deshpande and Montanari 
(2013) presented an iterative approximate message passing algorithm to solve the latter problem with 
sharp boundaries on SNR. However, in contrast to submatrix localization, where the signal submatrix 
can be located anywhere within the matrix, the clique Ending problem requires the signal to be centered 
on the diagonal. 

We would like to emphasize the difference between detection and localization problems. When M is 
a vector, Donoho and Jin (2004) proposed the “higher criticism” approach to solve the detection problem 
under the Gaussian sequence model. Combining the results in (Donoho and Jin, 2004; Ma and Wu, 
2013), in the computationally efficient region, there is no loss in treating M in model (2) as a vector and 
applying the higher criticism method to the vectorized matrix for the problem of submatrix detection, 
fn fact, the procedure achieves sharper constants in the Gaussian setting. However, in contrast to the 
detection probiem, we will show that for localization, it is crucial to utilize the matrix structure, even in 
the computationally efficient region. 
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1.3 Notation 

Let [m\ denote the index set {1,2,...,m}. For a matrix X £ U mxn , Xi. £ denotes its z'-th row and 
X.j e IR m denotes its y-th column. For any / £ [m],J Q \n], Xjj denotes the submatrix correspond¬ 
ing to the index set / x /. For a vector v £ IR", \\v\\( p = (Lze[n] \vi\ p ) 1,p and for a matrix M £ R mxn , 

||M ||£ = sup^^Q \\Mv\\e p l\\ v\\e . When p - 2 , the latter is the usual spectral norm, abbreviated as ||M|| 2 . 
The nuclear norm a matrix M is defined as a convex surrogate for the rank, with the notation to be || M|| *. 
The Frobenius norm of a matrix M is defined as \\M\\f = M?^.. The inner product associated with 

the Frobenius norm is defined as < A,B ) = tr(A r f>). 

Denote the asymptotic notation a{ri) = 0(fo(zz)) if there exist two universal constants c/, c u such that 
ci< lim a(n)/ b(n) < lim a(n)/b(n) < c„. 0* is asymptotic equivalence hiding logarithmic factors in the 

n—> 00 n^oo 

following sense: a{n) = 0*(h(n)) iff there exists c > 0 such that a[n) = 0(h(n)log c n). Additionally, we 
use the notation a{n) ~ b{ri) as equivalent to a{ri) = 0 (h(n)), a{ri) ^ b{n) iff lim^^oo a{ri)lb(ri) = 00 and 
a{n ) ^ b{n) iff lim^oo a{ri) / b{n) = 0 . 

We define the zero-mean sub-Gaussian random variable z with sub-Gaussian parameter o in terms 
of its Laplacian. If there exists a universal constant c > 0, 

te Az < exp(cr 2 A 2 /2c), forallA>0 (4) 

then we have 

P(|z| > erf) < 2 -exp(-c- 1 2 12). 

We call a random vector Z £U n isotropic with parameter a if 

Uv T Z) 2 = o 2 \\v\\) 2 , for all v£U n . 

Clearly, Gaussian and Bernoulli measures, and more general product measures of zero-mean sub-Gaussian 
random variables satisfy this isotropic definition up to a constant scalar factor. 


1.4 Our Contributions 


To state our main results, let us first define a hierarchy of algorithms in terms of their worst-case running 
time on instances of the submatrix localization problem: 

LinAIg c PolyAIg c ExpoAIg c AIIAIg. 

The set LinAIg contains algorithms sd that produce an answer (in our case, the localization subset Cff) 
in time linear in m x n (the minimal computation required to read the matrix). The classes PolyAIg and 
ExpoAIg of algorithms, respectively, terminate in polynomial and exponential time, while AIIAIg has no 
restriction. 

Combining Theorem 3 and 4 in Section 2 and Theorem 5 in Section 3, the statistical and computa¬ 
tional boundaries for submatrix localization can be summarized as follows. 


Theorem 1 (Computational and Statistical Boundaries). Consider the submatrix localization problem 
under the model ( 2 ). The computational boundarySNR c for the dense case whenmm{k m , k n } £3 max{ m 112 , n U2 \ 
is 


SNR r 


' m v n 
km kn 


■ + ■ 


I log n ^ log m 


(5) 
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in the sense that 


A 

lim inf sup P [r^ t R rn or C„ ^ C„) = 0 , if - >3 SNR c ( 6 ) 

m,n,k m ,k n —> oo^ e LinAlg ^ > O 

A 

lim inf sup P ^ R m or C„ # CJ > 0, i/'-^SNR c (7) 

m,n,k m ,k n -~ oo^ePolyAIg MeJ{ tJ 

where (7) holds under the Hidden Clique hypothesis HQ (see Section 2.1). For the sparse case when 
ma x{k m , k n } Q min {m 112 , n 112 }, the computational boundary zsSNR c = 0 * ( 1 ), more precisely 


1^SNR C ^ 



m v n 
km kn 


The statistical boundary SNR s is 


SNR s - 



logm 



in the sense that 


lim inf sup P( ^ R m orC„ ^ C„] = 0, 

m,7i,k m ,k n —^ °°,e/eExpoAlg MeJZ ^ 

lim inf sup P (r% ^ R m orC„ t C„] > 0, 

m,n,k m ,k n —>oo&dEM\Mg 


( 8 ) 


if ~(Z SNR s (9) 

a 

if — ^ SNR s (10) 

a 


under the minimal assumption ma x{k m , k n ] ^ min {m, n}. 

If we parametrize the submatrix model as m = n,k m — k n ~ k = 0*(n“),A/ff = 0*(n _ Q, for some 
0 < a, f < 1, we can summarize the results of Theorem 1 in a phase diagram, as illustrated in Figure 1. 



Figure 1: Phase diagram for submatrix localization. Red region (C): statistically impossible, where even 
without computational budget, the problem is hard. Blue region (B): statistically possible but compu¬ 
tationally expensive (under the hidden clique hypothesis), where the problem is hard to all polynomial 
time algorithm but easy with exponential time algorithm. Green region (A): statistically possible and 
computationally easy, where a fast polynomial time algorithm will solve the problem. 
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To explain the diagram, consider the following cases. First, the statistical boundary is 


llogn ^ log m 


which gives the line separating the red and the blue regions. For the dense regime a > 1/2, the computa¬ 
tional boundary given by Theorem 1 is 


/ m v n / log n log m 

y kmkn V k m k n 

which corresponds to the line separating the blue and the green regions. For the sparse regime a < 1/2, 
the computational boundary is 0(1) ^ SNR C ^ 0(y ; log p^-), which is the horizontal line connecting 
(a = 0,/3 = 0 )to(a = l/2,/3 = 0). 

As a key part of Theorem 1, we provide various linear time spectral algorithms that will succeed in 
localizing the submatrix with high probability in the regime above the computational threshold. Fur¬ 
thermore, the method is adaptive: it does not require the prior knowledge of the size of the submatrix. 
This should be contrasted with the method of Chen and Xu (2014) which requires the prior knowledge 
of k m , lc n ; furthermore, the running time of their SDP-based method is superlinear in nm. Under the 
hidden clique hypothesis, we prove that below the computational threshold there is no polynomial time 
algorithm that can succeed in localizing the submatrix. This is a new result that has not been established 
in the literature. We remark that the computational lower bound for localization requires a technique 
different from the lower bound for detection ; the latter has been resolved in Ma and Wu (2013). 

Beyond localization of one single submatrix, we generalize both the computational and statistical 
story to a growing number of submatrices in Section 2.5. As mentioned earlier, the statistical bound¬ 
ary for one single submatrix localization has been investigated by Butucea et al. (2013) in the Gaussian 
case. Our result focuses on the computational intrinsic difficulty of localization for a growing number of 
submatrices, at the expense of not providing the exact constants for the thresholds. 

The phase transition diagram in Figure 1 for localization should be contrasted with the correspond¬ 
ing result for detection, as shown in (Butucea and Ingster, 2013; Ma and Wu, 2013). For a large enough 
submatrix size (as quantified by a > 2/3), the computationally-intractable-but-statistically-possible re¬ 
gion collapses for the detection problem, but not for localization. In plain words, detecting the presence 
of a large submatrix becomes both computationally and statistically easy beyond a certain size, while 
for localization there is always a gap between statistically possible and computationally feasible regions. 
This phenomenon also appears to be distinct to that of other problems like estimation of sparse princi¬ 
pal components (Cai et al., 2013), where computational and statistical easiness coincide with each other 
over a large region of the parameter spaces. 


1.5 Organization of the Paper 

The paper is organized as follows. Section 2 establishes the computational boundary, with the compu¬ 
tational lower bounds given in Section 2.1 and upper bound results in Sections 2.2-2.4. An extension 
to the case of multiple submatrices is presented in Section 2.5. The upper and lower bounds for statis¬ 
tical boundary for multiple submatrices are discussed in Section 3. A discussion is given in Section 4. 
Technical proofs are deferred to Section 5. In addition to the spectral method given in Section 2.2 and 
2.4, Appendix A contains a new analysis of a known method that is based on a convex relaxation (Chen 
and Xu, 2014). Comparison of computational lower bounds for localization and detection is included in 
Appendix B. 
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2 Computational Boundary 


We characterize in this section the computational boundaries for the submatrix localization problem. 
Sections 2.1 and 2.2 consider respectively the computational lower bound and upper bound. The com¬ 
putational lower bound given in Theorem 2 is based on the hidden clique hypothesis. 

2.1 Algorithmic Reduction and Computational Lower Bound 

Theoretical Computer Science identifies a range problems which are believed to be “hard,” in the sense 
that in the worst-case the required computation grows exponentially with the size of the problem. Faced 
with a new computational problem, one might try to reduce any of the “hard” problems to the new 
problem, and therefore claim that the new problem is as hard as the rest in this family. Since statistical 
procedures typically deal with a random (rather than worst-case) input, it is natural to seek token prob¬ 
lems that are believed to be computationally difficult on average with respect to some distribution on 
instances. The hidden clique problem is one such example (for recent results on this problem, see Feld¬ 
man et al. (2013); Deshpande and Montanari (2013)). While there exists a quasi-polynomial algorithm, 
no polynomial-time method (for the appropriate regime, described below) is known. Following several 
other works on reductions for statistical problems, we work under the hypothesis that no polynomial¬ 
time method exists. 

Let us make the discussion more precise. Consider the hidden clique model c S{N,k) where N is the 
total number of nodes and k is the number of clique nodes. In the hidden clique model, a random 
graph instance is generated in the following way. Choose k clique nodes uniformly at random from all 
the possible choices, and connect all the edges within the clique. For all the other edges, connect with 
probability 1/2. 

Hidden Clique Hypothesis for Localization (H C|) Consider the random instance of hidden clique model 
^{N,k). For any sequence k(N) such that k[N) < NP for some 0 < yS < 1/2, there is no randomized poly¬ 
nomial time algorithm that can find the planted clique with probability tending to 1 as IV —> oo. Mathe¬ 
matically, define the randomized polynomial time algorithm class PolyAIg as the class of algorithms .<A 
that satisfies 


lim sup EciiqueP^uv.jOiciique (runtime of si not polynomial in N) = 0. 

JV,K(iV)^oo^ e p 0 | yA | g 


Then 

lim inf EciiquelP^(iv,K)|Ciique (clique set returned by sd not correct) > 0, 

JV,K(AO^oo^£PolyAlg 

where P > ^pv,K)|Ciique is the (possibly more detailed due to randomness of algorithm) cr-field conditioned 
on the clique location and Eciique is with respect to uniform distribution over all possible clique locations. 

Hidden Clique Hypothesis for Detection (HC c |) Consider the hidden clique model ( ${N,k). For any 
sequence of k(N) such that k(N) < for some 0 < /3 < 1/2, there is no randomized polynomial time 
algorithm that can distinguish between 


Ho : 2? er v.s. Y\ a : &> H c 
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with probability going to 1 as N — oo. Here , c 3%r is the Erdos-Renyi model, while S^hc is the hidden 
clique model with uniform distribution on all the possible locations of the clique. More precisely, 


lim inf EciiqueP’^(iv,K)|Ciique (detection decision returned by sd wrong) > 0 , 

JV,k(A0^oo,s/e PolyAIg 

where P^(jv,K)|Ciique and Eciique are the same as defined in HQ. 

The hidden clique hypothesis has been used recently by several authors to claim computational in¬ 
tractability of certain statistical problems. In particular, Berthet and Rigollet (2013a); Ma and Wu (2013) 
assumed the hypothesis HQ and Wang et al. (2014) used HQ. Localization is harder than detection, in 
the sense that if an algorithm st solves the localization problem with high probability, it also correctly 
solves the detection problem. Assuming that no polynomial time algorithm can solve the detection prob¬ 
lem implies impossibility results in localization as well. In plain language, H Q is a milder hypothesis than 
HC d . 

We will provide two computational lower bound results, one for localization and the other for de¬ 
tection, in Theorems 2 and 6 . The latter one will be deferred to Appendix B to contrast the difference of 
constructions between localization and detection. The detection computational lower bound was first 
proved in Ma and Wu (2013). For the localization computational lower bound, to the best of our knowl¬ 
edge, there is no proof in the literature. Theorem 2 ensures the upper bound in Lemma 1 being sharp. 

Theorem 2 (Computational Lower Bound for Localization). Consider the submatrix model (2) with pa¬ 
rameter tuple [m = n,k m - k n - n“,A/cr = n~h, where ~ < a < 1, f > 0. Under the computational as¬ 
sumption HQ, if 



knxkyi 


p > a 


1 

2 ’ 


it is not possible to localize the true support of the submatrix with probability going to 1 within polynomial 
time. 


Our algorithmic reduction for localization relies on a bootstrapping idea based on the matrix struc¬ 
ture and a cleaning-up procedure introduced in Lemma 13 given in Section 5. These two key ideas offer 
new insights in addition to the usual computational lower bound arguments. Bootstrapping introduces 
an additional randomness on top of the randomness in the hidden clique. Careful examination of these 
two cr-fields allows us to write the resulting object into mixture of submatrix models. For submatrix 
localization we need to transform back the submatrix support to the original hidden clique support ex¬ 
actly, with high probability. In plain language, even though we lose track of the exact location of the 
support when reducing the hidden clique to submatrix model, we can still recover the exact location of 
the hidden clique with high probability. For technical details of the proof, please refer to Section 5. 


2.2 Adaptive Spectral Algorithm and Computational Upper Bound 

In this section, we introduce linear time algorithm that solves the submatrix localization problem above 
the computational boundary SNR C . Our proposed localization Algorithms 1 and 2 is motivated by the 
spectral algorithm in random graphs (McSherry, 2001; Ng et al., 2002). 
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Algorithm 1: Vanilla Spectral Projection Algorithm for Dense Regime 
Input: X e M mxn the data matrix. 

Output: A subset of the row indexes R m and a subset of column indexes C n as the localization sets 
of the submatrix. 

1. Compute top left and top right singular vectors U.\ and Vfi, respectively (these correspond to 
the SVD X = U1V T ); 

2. To compute C n , calculate the inner products 1 < j < n. These values form two clusters. 

Similarly, for the R m , calculate X,-. V.\, 1 < i < m and obtain two separated clusters. A simple 
thresholding procedure returns the subsets C n and R m . 


The proposed algorithm has several advantages over the localization algorithms that appeared in 
literature. First, it is a linear time algorithm (that is, 0(mn) time complexity). The top singular vectors 
can be evaluated using fast iterative power methods, which is efficient both in terms of space and time. 
Secondly, this algorithm does not require the prior knowledge of k m and k„ and automatically adapts to 
the true submatrix size. 

Lemma 1 below justifies the effectiveness of the spectral algorithm. 

Lemma 1 (Guarantee for Spectral Algorithm). Consider the submatrix model (2), Algorithm 1 and assume 
min{fc OT , k n } X max{m 1/2 , n 112 }. There exist a universal C > 0 such that when 


>C- 


o 


imy n 
km kn 


■ + ■ 


llogn log m 


the spectral method succeeds in the sense that R m = R m , C n = C n with probability at least 1 - m c - n c - 
2exp(-c(m + n)). 


2.3 Dense Regime 

We are now ready to state the SNR boundary for polynomial-time algorithms (under an appropriate 
computational assumption), thus excluding the exhaustive search procedure. The results hold under 
the dense regime when k X n 112 . 

Theorem 3 (Computational Boundary for Dense Regime). Consider the submatrix model (2) and assume 
min{k OT , k n } X max{m 1/2 , n 112 }. There exists a critical rate 


SNR C = 


imy n 
km k tl 


■ + 



logm 



for the signal to noise ratio SNR c such that for XI o "f SNR C , both the adaptive linear time Algorithm 1 and 
the robust polynomial time Algorithm 5 will succeed in submatrix localization, i.e., R m = R m ,C n = C n , 
with high probability. For XIa f SNR C , there is no polynomial time algorithm that will work under the 
hidden clique hypothesis HC|. 

The proof of the above theorem is based on the theoretical justification of the spectral Algorithm 1 
and convex relaxation Algorithm 5, and the new computational lower bound result for localization in 
Theorem 2. We remark that the analyses can be extended to multiple, even growing number of subma¬ 
trices case. We postpone a proof of this fact to Section 2.5 for simplicity and focus on the case of a single 
submatrix. 


9 

















2.4 Sparse Regime 

Under the sparse regime when k ^ n 112 , a naive plug-in of Lemma 1 requires the SNR C to be larger 
than ®(n ll2 lk) £3 \/\ogn, which implies the vanilla spectral Algorithm 1 is outperformed by simple en- 
trywise thresholding. However, a modified version with entrywise soft-thresholding as a preprocessing 
de-noising step turns out to provide near optimal performance in the sparse regime. Before we introduce 
the formal algorithm, let us define the soft-thresholding function at level t to be 


t 7 r(y) = sign(y)(|y|-f) + . ( 11 ) 

Soft-thresholding as a de-noising step achieving optimal bias-and-variance trade-off has been widely 
understood in the wavelet literature, for example, see Donoho and Johnstone (1998). 

Now we are ready to state the following de-noised spectral Algorithm 2 to localize the submatrix 
under the sparse regime when k ^ n 1/2 . 


Algorithm 2: De-noised Spectral Algorithm for Sparse Regime 
Input: X e M' nxn the data matrix, a thresholding level t = 0( a y'log jfff). 

Output: A subset of the row indexes R m and a subset of column indexes C„ as the localization sets 
of the submatrix. 

1. Soft-threshold each entry of the matrix X at level t, denote the resulting matrix as f] fX) ; 

2. Compute top left and top right singular vectors U.\ and Hi of matrix r] t (X), respectively (these 
correspond to the SVD t] t [X) = UZV T ); 

3. To compute C n , calculate the inner products U T X ■ t]r(X.j), 1 < j <n. These values form two 
clusters. Similarly, for the R m , calculate r/ f (A ; -.) • V.\, 1 < i < m and obtain two separated clusters. A 
simple thresholding procedure returns the subsets C n and R m . 


Lemma 2 below provides the theoretical guarantee for the above algorithm when k^n 112 . 


Lemma2 (Guarantee for De-noised Spectral Algorithm). Consider the submatrix model (2), soft-thresholded 
spectral Algorithm 2 with thresholded levelot, and assume min{L,„, k n } 33 max{m 1/2 , n U2 \. There exist a 
universal C > 0 such that when 


>C- 


m v n /log n log m 

k m kn + VXr V- fcT 


■e~ t2,2 + t 


the spectral method succeeds in the sense thatR m = R m , C n = C n with probability at least 1 - m~ c - n~ c - 
2exp(-c(m +n)). Further if we choose t = ®[a ogf^p) as the optimal thresholding level, we have de- 
noised spectral algorithm works when 

A / mv n 

a-]/ l 0g kJ^- 

Combining the hidden clique hypothesis H C| together with Lemma 2 , we have the following theorem 
holds under the sparse regime when k 33 n 112 . 


Theorem 4 (Computational Boundary for Sparse Regime). Consider the submatrix model (2) and assume 
ma x{k m , k n \ f min {m 1/2 , n 112 }. There exists a critical rate for the signal to noise ratio SN R c between 


1^SNR C ^ 



m v n 
km kn 
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such that for XIa £3 \/lof> jpF"> the linear time Algorithm 2 will succeed in submatrix localization, i.e., 

R m = R m ,C n = C n , with high probability. For Ala 33 1 , there is no polynomial time algorithm that will 
work under the hidden clique hypothesis HC|. 

Remark 4.1. The upper bound achieved by the de-noised spectral Algorithm 2 is optimal in the two 
boundary cases: k = 1 and k ~ n 112 . When k = 1, both the information theoretic and computational 
boundary meet at ylog n. When k - n 112 , the computational lower bound and upper bound match in 
Theorem 4, thus suggesting the near optimality of Algorithm 2 within the polynomial time algorithm 
class. The potential logarithmic gap is due to the crudeness of the hidden clique hypothesis. Precisely, 
for k = 2 , hidden clique is not only hard for G(n, p) with p = 1 / 2 , but also hard for G(», p) with p = 1 /logn. 
Similarly for k = n a , a < 1/2, hidden clique is not only hard for G(n, p) with p = 1/2, but also for some 

0<p< 1/2. 


2.5 Extension to Growing Number of Submatrices 

The computational boundaries established in the previous sections for a single submatrix can be ex¬ 
tended to non-overlapping multiple submatrices model (3). The non-overlapping assumption corre¬ 
sponds to that for any l<s^t<r,R s nR t = 0 and C s n C f = 0. The Algorithm 3 below is an extension of 
the spectral projection Algorithm 1 to address the multiple submatrices localization problem. 

Algorithm 3: Spectral Algorithm for Multiple Submatrices 
Input: X £ [R ,nx “ the data matrix. A pre-specified number of submatrices r. 

Output: A subset of the row indexes {R s m , 1 < 5 < r] and a subset of column indexes {C s n , 1 < s < r} 
as the localization of the submatrices. 

1. Calculate top r left and right singular vectors in the SVD X = U2.V 1 . Denote these vectors as 
U r £ U mxr and V r £ U nxr , respectively; 

2. For the C s n , 1 < s < r, calculate the projection U r [UjU r )~ 1 Uf X.j, 1 < j < n, run /c-means 
clustering algorithm (with k = r + 1) for these n vectors in R" 1 . For the R s m , 1 < s<r, calculate 
V r {Vj V r )~ l VjxT, \<i < m, run /c-means clustering algorithm (with k = r + 1) for these m 
vectors in IR“ (while the effective dimension is IR r ). 


We emphasize that the following Proposition 3 holds even when the number of submatrices r grows 
with m, n. 


Lemma 3 (Spectral Algorithm for Non-overlapping Submatrices Case). Consider the non-overlapping 
multiple submatrices model (3) and Algorithm 3. Assume 

k^-k m ,k^~k n ,\ s ~X 

for all 1 < s< r and min{/c m , k n } Q max{m 1/2 , n 112 }. There exist a universal C > 0 such that when 


A 

->C- 

cr 


km A k n 


■ + ' 


I log n / log m 


kn 


kn 


+ ' 


1 m v n 


km k n 


( 12 ) 


the spectral method succeeds in the sense that R^ = , C, ( , s) = C )/ 1 ,1 < s < r with probability at least 

1 - m~ c - n~ c -2 exp (-c(m + n)). 


Remark 4.2. Under the non-overlapping assumption, rk m f m, rk n f n hold in most cases. Thus the 
first term in Equation (12) is dominated by the latter two terms. Thus a growing number r does not affect 
the bound in Equation (12) as long as the non-overlapping assumption holds. 
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3 Statistical Boundary 


In this section we study the statistical boundary. As mentioned in the introduction, in the Gaussian noise 
setting, the statistical boundary for a single submatrix localization has been established in Butucea et al. 
(2013). In this section, we generalize to localization of a growing number of submatrices, as well as sub- 
Gaussian noise, at the expense of having non-exact constants for the threshold. 


3.1 Information Theoretic Bound 

We begin with the information theoretic lower bound for the localization accuracy. 

Lemma 4 (Information Theoretic Lower Bound). Consider the submatrix model (2) with Gaussian noise 
Zjj ~ -jY({),(J 2 ). For any fixed 0 < a < 1, there exist a universal constant C a such that if 


A ^ llogiml k m ) log {nlk n ) 

a Ca '\ - k - + - k -’ (13) 

u v K 'm 

any algorithm sd will fail to localize the submatrix with probability at least 1 -a - k i 0 g(»/t- j 

in the following minimax sense: 


inf sup 

■s/eANAIg MeM 



or C„ d- C n ] > 1 -a- 


_log2_ 

k m log {ml km) + k n log [n! k n )' 


3.2 Combinatorial Search for Growing Number of Submatrices 

Combinatorial search over all submatrices of size k m x k n finds the location with the strongest aggregate 
signal and is statistically optimal (Butucea et al., 2013; Butucea and Ingster, 2013). Unfortunately, it re¬ 
quires computational complexity 0 + ( k ) j, which is exponential in k m , k n . The search Algorithm 4 
was introduced and analyzed under the Gaussian setting for a single submatrix in Butucea and Ingster 
(2013), which can be used iteratively to solve multiple submatrices localization. 

Algorithm 4: Combinatorial Search Algorithm 
Input: X e K mxn the data matrix. 

Output: A subset of the row indexes R m and a subset of column indexes C„ as the localization of 
the submatrix. 

For all index subsets / x / with |/| = k m and |/| = k n , calculate the sum of the entries in the 
submatrix Xjj. Report the index subset R m x C n with the largest sum. 

For the case of multiple submatrices, the submatrices can be extracted with the largest sum in a 
greedy fashion. 

Lemma 5 below provides a theoretical guarantee for Algorithm 4 to achieve the information theoretic 
lower bound. 

Lemma 5 (Guarantee for Search Algorithm). Consider the non-overlapping multiple submatrices model 
(3) and iterative application of Algorithm 4 in a greedy fashion for r times. Assume 

k™~k m ,k™-k n , 
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for alll< s< r and rnax{/c m , k n ] ^ min {m, n}. There exists a universal constant C > 0 such that if 



then Algorithm 4 will succeed in returning the correct location of the submatrix with probability at least 


mn 


To complete Theorem 1, we include the following Theorem 5 capturing the statistical boundary. It is 
proved by exhibiting the information-theoretic lower bound Lemma 4 and analyzing Algorithm 4. 

Theorem 5 (Statistical Boundary). Consider the submatrix model (2). There exists a critical rate 



for the signal to noise ratio, such that for any problem with A / a £3 S N R s , the statistical search Algorithm 4 
will succeed in submatrix localization, i.e., R m = R m , C n = C n , with high probability. On the other hand, if 
Ala f SNR s , no algorithm will work (in the minimax sense) with probability tending to 1 . 

4 Discussion 

In this paper we established the computational and statistical boundaries for submatrix localization in 
the setting of a growing number of submatrices with subgaussian noise. The primary goals are to demon¬ 
strate the intrinsic gap between what is statistical possible and what is computationally feasible and to 
contrast the interplay between computational efficiency and statistical accuracy for localization with 
that for detection. 

Submatrix Localization v.s. Detection As pointed out in Section 1.4, for any k = n a ,0 < a < 1, there is 
an intrinsic SNR gap between computational and statistical boundaries for submatrix localization. Un¬ 
like the submatrix detection problem where for the regime 2/3 < a < 1, there is no gap between what 
is computationally possible and what is statistical possible. The inevitable gap in submatrix localiza¬ 
tion is due to the combinatorial structure of the problem. This phenomenon is also seen in some net¬ 
work related problems, for instance, stochastic block models with a growing number of communities. 
Compared to the submatrix detection problem, the algorithm to solve the localization problem is more 
complicated and the techniques required for the analysis are much more involved. 

Detection for Growing Number of Submatrices The current paper solves localization of a growing 
number of submatrices. In comparison, for detection, the only known results are for the case of a single 
submatrix as considered in Butucea and Ingster (2013) for the statistical boundary and in Ma and Wu 
(2013) for the computational boundary. The detection problem in the setting of a growing number of 
submatrices is of significant interest. In particular, it is interesting to understand the computational and 
statistical trade-offs in such a setting. This will need further investigation. 
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Estimation of the Noise Level a Although Algorithms 1 and 3 do not require the noise level a as an 
input, Algorithm 2 does require the knowledge of a. The noise level o can be estimated robustly. In 
the Gaussian case, a simple robust estimator of a is the following median absolute deviation (MAD) 
estimator due to the fact that M is sparse: 


<7 = median,- j |- median,(X,y) | /O 1 (0.75) 
~ 1.4826 x median|- median,y (X,y)|. 


5 Proofs 

We prove in this section the main results given in the paper. We first collect and prove a few important 
technical lemmas that will be used in the proofs of the main results. 

5.1 Prerequisite Lemmas 

We start with two Lemmas 6 and 7 that due to perturbation theory. 

Lemma 6 (Stewart and Sun (1990) Theorem 4.1). Suppose that A = A + E, all of which are matrices of the 
same size, and we have the following singular value decomposition 


Ii 0 


[U 1 ,U 2 ,U 3 ] T A[V l ,V 2 ] = 0 Z 2 

0 0 


(14) 


II 0 


[U 1 ,U 2 ,U3] T A[V 1 ,V 2 \= 0 z 2 

0 0 


(15) 


Let O be the matrix of canonical angles between S&{U\) and £&{U\), and let © be the matrix of canonical 
angles between £&{Vi) andS%(V\) (here LA denotes the linear space). Define 


R = AV\ — UiDi 
S = A t U 1 -V i£i 


(16) 

(17) 


Then suppose there is a numbers > 0 such that 


min|cr(Zi)-o r (Z 2 )| > S and mincr(Zi) > <5. 


Then 



Further, suppose there are numbers a, 8 such that 


mincr(Zi) > 5 + a and maxcr(Z 2 ) < a, 


then for 2-norm, or any unitarily invariant norm, we have 


max{|| sin<f)|| 2 , || sin©|| 2 }< 


max{||J?|| 2 ,||S|| 2 } 

8 
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Let us use the above version of the perturbation bound to derive a lemma that is particularly useful 
in our case. Simple algebra tells us that 


A[V l ,V 2 ] = \U 1 ,U2,U 3 ] 


±i 0 
0 1 2 
0 0 


= [Uili,U 2 l 2 ]. 


(18) 


A[V 1i V 2 ] = [AV 1 ,AV 2 ] (19) 

LA-A) [Vi, V 2 ] = [Uil i - AV x ,U 2 ±2 - AV 2 ] (20) 

\\A-A\\ 2 F = Tr{lA-A)[V 1 ,V 2 ][V 1 ,V 2 ] T {A-A) T ) (21) 

= \\R\\ 2 p + \\U 2 1 2 -AV 2 \\ 2 f >\\R\\ 2 p . (22) 

Similarly, we have 

[A- A) T [U\, [j 2> U 3 \ = [Vi£i - A r U lt V 2 1 2 - A t U 2 ,-A t U 3 \ (23) 

\\A-A\\ 2 F = Tr((A-A) T [U 1 ,U 2 ,U 3 \[U l ,U 2 ,U 3 ] T (A-A)) (24) 

= || S||| + || V 2 Z 2 - A t U 2 ||| + || A t U 3 ||| > || S|||. (25) 


Thus, it holds that 


\\A - A || F > max( || R || F , || S || F ) 


and similarly we have (since the operator norm of a whole matrix is larger than that of the submatrix) 


||A-A|| 2 >max(||i?|| 2 ,||S|| 2 ). 


Thus the following version of the Wedin’s Theorem holds. 

Lemma 7 (Davis-Kahan-Wedin’s Type Perturbation Bound). It holds that 


sin<I>||! + || sin 0 ||| < 


V2\\E\\j 


" F " -" F “ 5 

and also the following holds for 2-norm (or any unitary invariant norm) 

\\E\\ 2 

max{|| sinO|| 2 , || sin 0 || 2 } < ——. 


We will then introduce some concentration inequalities. Lemmas 8 and 9 are concentration of mea¬ 
sure results from random matrix theory. 


Lemma 8 (Vershynin (2010), Theorem 39). Let Z e u mxn be a matrix whose rows Z,-. are independent 
sub-Gaussian isotropic random vectors in R" with parameter a. Then for every t> 0, with probability at 
leastl-2exp{-ct 2 ) one has 

\\Z\\ 2 < a{\fm + C\fn+ t ) 
where C, c > 0 are some universal constants. 
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Lemma 9 (Hsu et al. (2012), Projection Lemma). Assume Z e R n is an isotropic sub-Gaussian vector with 
i.i.d. entries and parameter a. PA is a projection operator to a subspace of dimension r, then we have the 
following concentration inequality 

PGI^ZH^ > a 2 {r + 2\/7t + 2t)) < exp(-cf), 
where c> 0 is a universal constant. 


The proof of this lemma is a simple application of Theorem 2.1 inHsu et al. (2012) for the case that SA 
is a rank r positive semidefinite projection matrix. 

The following two are standard Chernoff-type bounds for bounded random variables. 


Lemma 10 (Hoeffding (1963), Hoeffding’s Inequality). LetXi, 1 < i < nbe independent random variables. 
Assume a; < < bi,\ < i < n. Then for S n = £” =1 X{ 


P(|S n 


ES n \ > u) < 2 exp 


2 IT 


ru ibr 


(26) 


Lemma 11 (Bennett (1962), Bernstein’s Inequality). Let X/, 1 < i < n be independent zero-mean random 
variables. Suppose \X- L \ < M, 1 < i < n. Then 


1 n ' 


\Y J X i> u 

V«=t ) 

< exp 


u 2 / 2 

I'L, EX ; 2 + Mu/3 


(27) 


We will end this section stating the Fano’s information inequality, which plays a key role in many 
information theoretic lower bounds. 


Lemma 12 (Tsybakov (2009) Corollary 2.6). LetfAo,fA],...,fA^ be probability measures on the same prob¬ 
ability space (0,^), M > 2. If for some 0 < a < 1 


1 

M+ 1 


M 


X rf KL i&i 11 ^ a ■ log M 
i =0 


where 


Then 


M 


SA = 


M + 


1 i =0 


log(M+ l)-log 2 

p e M — Pe M — - ;- Ct 

y ’ ' log M 


where p Bi M is the minimax error for the multiple testing problem. 


(28) 


(29) 


5.2 Main Proofs 

Proof of Lemma 1. Recall the matrix form of the submatrix model, with the SVD decomposition of the 
mean signal matrix M 

X = A y/k m k n UV T + Z. 

The largest singular value of A UV T is A\/k m lc n , and all the other singular values are 0s. Davis-Kahan- 
Wedin’s perturbation bound tells us how close the singular space of X is to the singular space of M. 


16 













Let us apply the derived Lemma 7 to X = A\Jk m k n UV T + Z. Denote the top left and right singular 
vector of X as U and V. One can see that E||Z ||2 - cr(\/m + yfn) under very mild hnite fourth moment 
conditions through a result in (Latala, 2005) . Lemma 8 provides a more explicit probabilisitic bound for 
the concentration of the largest singular value of i.i.d sub-Gaussian random matrix. Because the rows 
Z/. are sampled from product measure of mean zero sub-Gaussians, they naturally satisfy the isotropic 
condition. Hence, with probability at least 1 - 2exp (-c(m + ri)), via Lemma 8, we reach 

\\Z\\ 2 < C-cr{\/m + \/n). (30) 

Using Weyl’s interlacing inequality, we have 

\adX)-aim\<mz 


and thus 


0-1CX") > A\Jk m k n - II Z|| 2 
o- 2 (X)<||Z|| 2 . 


Applying Lemma 7, we have 

max{|sinZ([7, t/)|, | sinZ(V, V)|} < 

In addition 


Co{\fm+ \fn) 


_ o{\fm + \[ri) 

\\Jk m k n - Ca{s/m + sfn) A\Jk m k n 


U-Ulg, = \/2-2cosZ[U,U) = 2\sm-Z{U,U)\, 


which means 


max{||[/- U\\( 2> \\V - V\\f 2 } < C 


cr{s/m + \/n) 

A\J k m k n 


And according to the definition of the canonical angles, we have 


max{ || UU T - UU t \\ 2 , II VV T - VV T || 2 } < C 


cr( v / m + y/n) 
A \J k m k n 


Now let us assume we have two observations of X. We use the first observation X to solve for the 
singular vectors U, V, we use the second observation X to project to the singular vectors U, V. We 
can use Tsybakov’s sample cloning argument (Tsybakov (2014), Lemma 2.1) to create two independent 
observations of X when noise is Gaussian as follows. Create a pure Gaussian matrix Z' and define 
X\ = X + Z' - M + (Z + Z') and X 2 = X- Z' = M + (Z - Z 1 ), making X\,X 2 independent with the vari¬ 
ance being doubled. This step is not essential because we can perform random subsampling as in Vu 
(2014); having two observations instead of one does not change the picture statistically or computation¬ 
ally. Recall X = M+ Z = A ^/k m k„UV T + Z. 

Define the projection operator to be SZ, we start the analysis by decomposing 


\\9> t jX. j -M. } \\ e2 <\\& t j{X. j 


M.j) || + II ~ I) M.j || g 2 


(31) 


for 1 < j < n. 


17 














For the first term of (31), note that X.j - M.j = Z.j e IR' n is an i.i.d. isotropic sub-Gaussian vector, and 
thus we have through Lemma 9, for t = (1 + 1 / c) log n, Z.j e IR m , 1 < j < n and r = 1 


- M.j) \\( 2 >a\fr \ 1 + 2V1 + 1/C-- 
We invoke the union bound for all 1 < j < n to obtain 


/ log n , log n 

——+ 2(l + l/c)- —— 


< n 


-c -1 


(32) 


max WSZ^X.j - M.j) \\ e < o 

1 <j<n 1 1 


\[r + \/2(\ + \l c)-o\J\og 


n 


(33) 


< o + C- a 



(34) 


with probability at least 1 - n~ c . 

For the second term M.j = X.j - Z.j of (31), there are two ways of upper bounding it. The first ap¬ 
proach is to split 


|| - DM\\ 2 < || (S» & - I)X|| 2 + || - I)Z\\ 2 < 2||Z|| 2 . (35) 

The first term of (35) is cr 2 (X] < cr 2 {M) + ||Z ||2 through Weyl’s interlacing inequality, while the second 
term is bounded by ||Z|| 2 . We also know that ||Z|| 2 < C 2 -o{\/m + \/~n). Recall the definition of the induced 
P. 2 norm of a matrix (PZq - I)M: 

H2? 0 -I)M\\ 2 > —^ nxVKAuh, > VbtU&o-nM-M- 

II v lh 2 


In the second approach, the second term of (31) can be handled through perturbation Sin Theta Theorem 
7: 

•r t „,^|| _ ^ o\/m + n 


|| (£^p - 1)M.j || e2 = || - &u)M.j h 2 < || UU 1 - UU 1 || 2 • || M.j ||/ 2 < C- 


d."\/ k m . 


X \J k m k n 

This second approach will be used in the multiple submatrices analysis. 

Combining all the above, we have with probability at least 1 - n~ c - m~ c , for all 1 < j < n 


\\&>jyX.j - M.j\\g 2 < C■ ^°y log^ + °^ 
Similarly we have for all 1 < i < m, 


1 m v n 




WSZyXi. -Mj.. II ^ < C-\aJ\ogm + a 


i m v n 


(36) 


(37) 


Clearly we know that for i e R m and i' e [m] \R m 

ml - Ml\\g 2 = 

and for j e C n and j' e [n\\C n 

\\M.j - M.j>\\g 2 =A s/km 
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Thus if 


^n/^ot ^ 6 C- crJ\ogn + a 


lmy n 

k n 


A \fk~n > 6 C • a \/log m + a 


imy n 
km 


(38) 


(39) 


hold, then we have learned a metric d (a one dimensional line) such that on this line, data forms clusters 
in the sense that 

2 max | dj - df \ < min | d,- - d,/ 1 . 

i,i'^R m iE.R m ,i'e[m]\R m 

In this case, a simple cut-off clustering recovers the nodes exactly. 

In summary, if 


A> C-o 


( 


log n /log m m + n 


kn 


■ + 


k„ 


+1 


km kn 


the spectral algorithm succeeds with probability at least 

1 - m~ c - n~ c - 2 exp (-c(m + n)). 


□ 


Proof of Lemma 2. The proof of the validity of thresholded spectral algorithm at level o t is easy based on 
the proof of Lemma 1. Firstly we have the following decomposition 


d<Jt(X) — M + rhj t (Z) + B 

= ^Rm l C n +T l<Jt(Z) +B 


where B is the bias matrix satisfying 

Bjj — 0, if (i, j ) <£ B-m x CA 
\Bij\<2at, if {i,j)eR m xC„. 

Let us prove this fact. Clearly if (i,j) € R m x C n , Bij = 0. If (i,j) e R m x C n , we have 

\Bij\ = \ri at (X + Zij)- A-Tj at {Zij)\ 

^ |f?<rr(A + Zij) — (A+ Zjj )| + |(A+ Zjj) — A — Zif + \Zij — rj at {Zi y)| 

< 2a t 

where the last step uses \rjtiy) - y\ < t, for any y. Let us bound the variance of each thresholded entry 

dot (Zij), 


9 r 

[dotiZij)] 2 = 

Jo 

-r 

-JC 


2z-2P(j] at (Zij) > z)dz 


4zP{Za > z + at)dz 


-t] 

°° a I (z + 

4zexpf-c 


+ ot) 2 1 

2u 2 | 


dz 


< C-cr-exp(-- 
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for some universal constant C. Clearly after thresholding, r} at (Z) still have i.i.d entries, but the variance 
has been significantly reduced as t —■ oo. 

Via the perturbation analysis established in Proof of Lemma 1 


WrjotiZ) + B\\2< \\r]at(Z)\\2 + II5II 2 


S C-a,/^.exp(- T ) + 2 v^« 

as B only have k m k n non zero entries. Thus applying Lemma 7, we have 


max{|sinZ([7, L0|,|sinZ(V, V)|} < 


C-o\/m v n-exp(-^) + 2y/k m k n at 
X\Jk m k n - C■ a v/m v n-exp(-y) - 2\Jk m k n ot 
yJmM n-crexp(-y) + \]k m k n ot 
4 \/ kf)i k n 


As usual, we continue the analysis by decomposing (following the steps as in Lemma 7, but with an 
additional bias term B ) 


\\&> t jT](X.j)-M.j \\ e2 < \\& 0 [r 1 [X.j)-M.j )\\ e2 + U& t y-I)M.j \\ £2 

< w&uTjiz.jn^ + \\B.j\\f 2 + n&o-nM-M 

\Jm\i n-o exp(-y) + yjk m k n at 




T \/ kin k n 

for 1 < j < n. We know for j e C n and j' e [ri]\C n 

II M.j - M.j' |h 2 = A \fk~m 

Thus if 

s/mv n-crex p(-y)+ y k m k„ot 


A\/ k m , 


A \fk^i > 6C• |crexplog n + \/ r k m at + 


\/~kri 


(40) 


hold, then we have learned a metric d (a one dimensional line) such that on this line, data forms clusters. 
In this case, a simple cut-off clustering recovers the nodes exactly. 

In summary, if 


- >C- 
a 


m v n / log n log m 


km kn 


■ + 


kn 


kn 


■e f2/2 + t 


the thresholded spectral algorithm succeeds with probability at least 

1 - m~ c - n~ c - 2exp (-c(m+ n)). 


□ 

Proof of Theorem 2. Computational lower bound for localization (support recovery) is of different nature 
than the computational lower bound for detection (two point testing). The idea is to design a random¬ 
ized polynomial time algorithmic reduction to relate a an instance of hidden clique problem to our sub¬ 
matrix localization problem. The proof proceeds in the following way: we will construct a randomized 
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polynomial time transformation 3~ to map a random instance of ( ${N,k) to a random instance of our 
submatrix Jiim = n,k m - k n - k,Ala) (abbreviated as M{n, k,A/cr)). Then we will provide a quanti¬ 
tative computational lower bound by showing that if there is a polynomial time algorithm that pushes 
below the hypothesized computational boundary for localization in the submatrix model, there will be a 
polynomial time algorithm that solves hidden clique localization with high probability (a contradiction 
to HC|). 

Denote the randomized polynomial time transformation as 

3~ : <^(A/,k(AT)) — M[n, k = n a ,A/a = n~P). 


There are several stages for the construction of the algorithmic reduction. First we define a graph c S e iN, k (A/)) 
that is stochastically equivalent to the hidden clique graph k(A 0), but is easier for theoretical anal¬ 
ysis. c S e has the property: each node independently has the probability k(N)IN to be a clique node, and 
with the remaining probability a non-clique node. Using Bernstein’s inequality and the inequality (46) 
proved below, with probability at least 1 - 2 AT 1 the number of clique nodes K e in c S e 


r 

k 1 

k 


I 4logA/ 


< K <K 



• K 


(41) 


as long as k £ log N. 

Consider a hidden clique graph c S e {2N,2K[N )) with N = n and k[N) = k. Denote the set of clique 
nodes for c S e [2N ,2 k{N)) to be Cn,k • Represent the hidden clique graph using the symmetric adjacency 
matrix G e {-1, \} 2Nx2N t where G,y = 1 if i,j e Cn,k> otherwise with equal probability to be either -1 or 
1. As remarked before, with probability at least 1 - 2AG 1 , we have planted 2 k( 1 + o(l)) clique nodes in 
graph c S e with 2N nodes. Take out the upper-right submatrix of G, denote as Gjjr where U is the index 
set 1 < i < N and R is the index set N + 1 < j < 2N. Now Gjjr has independent entries. 

The construction of ST employs the Bootstrapping idea. Generate l 2 (with l ~ nJ’,0 < /) < 1 12) matri¬ 
ces through bootstrap subsampling as follows. Generate l - 1 independent index vectors e R", 1 < 
s < l, where each element rj/ ls> (/'), 1 < i < n is a random draw with replacement from the row indices \n\. 
Denote vector 0 (O) (i) = i, 1 < i < n as the original index set. Similarly, we can define independently the 
column index vectors 1 < t < l. We remark that these bootstrap samples can be generated in polyno¬ 
mial time U(Z 2 n 2 ). The transformation is a weighted average of l 2 matrices of size n x n generated based 
on the original adjacency matrix Gjjr- 


ST: 


Mjj ( Gjjr )^/( s )(i)0f f )(y), 1 ^ i,j S n. 

^ 0<s,t<l 


(42) 


Due to the bootstrapping property, the matrices 


( i)(pu ) c j) 


, indexed by 0 < s, t < l are in- 

1 <i,j<n 


dependent of each other. Recall that Cjv.k stands for the clique set of the hidden clique graph. We define 
the row candidate set Rj := {i e [n] : 3 0 < s < ky/^Hi) e Cjv.kI and column candidate set C/ := {j £ [n\ : 
3 0 < t < l, <p (r> ( j ) e C n , k }. Observe that Rj x C/ are the indices where the matrix M contains signal. 

There are two cases for M,y, given the candidate set Rj x C/. If i e R t and j e C/, namely when (/, j) is 
a clique edge in at least one of the l z matrices, then E[Mjy |^ e ] > l~ l where the expectation is taken over 
the bootstrap cr-field conditioned on the candidate set R[ x C/ and the original cr-field of c S e . Otherwise 
E[Mjj\ ( S e ] = K n ^} k 2 - for (i,j) t Ri* C/, where \E\ is a Binomial(Af 2 -k 2 , 1/2). With high probability, 
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E[M/j|^ e ] - -y= 2 ~ j t = o(j). Thus the mean separation between the signal position and non-signal 

position is \ ~ j 7 - Note in the submatrix model, it does not matter if the noise has mean zero or not 
(since we can subtract the mean)- only the signal separation matters. 

Now let us discuss the independence issue in M through our Bootstrapping construction. Clearly 
due to sampling with replacement and bootstrapping, condition on c S e , we have independence among 
samples for the same location (i,j) 




For the independence among entries in one Bootstrapped matrix, clearly 


-L (GuR)yW{i<)<f,W{j'y 

The only case where there might be a slight dependence is between {Gur) y/^u)(l> u Hj) and {Gur )^^. 
The way to eliminate the slight dependence is through Vu (2008) ’s result on universality of random dis¬ 
crete graphs. Vu (2008) showed random regular graph ‘Sin, n/2) shares many similarities as Erdos-Renyi 
random graph ‘Sin, 1 / 2 ), for instance, top and second eigenvalues (n /2 and \fn respectively), limiting 
spectral distribution, sandwich conjecture, determinant, etc. Let us consider the case where the upper- 
right of the adjacency matrix G consists of random bi-regular graph (see Deshpande and Montanari 
(2013) for difficulty of clique problem under random regular graph) with degree n/2 instead of the Erdos- 
Renyi graph. The only thing we need to change is assuming hidden clique hypothesis is still valid for the 
following random graph: for a n x n adjacency matrix G, first find a clique/principal submatrix of size k 
uniformly randomly and connect density, for the remaining part of the matrix, sample a random regular 
graph of G(n - k, ^jp) and a random bi-regular graph of size k x (n - k) with left regular degree n/2 - k 
and right regular degree Id2 (here degree test will not work in this graph and spectral barrier still suggests 
k d yfn is hard due to universality result of random discrete graphs). In the bootstrapping step, condi¬ 
tion on the same row i/c w (z) being not a clique, ( Gur -L (G[/i?)wW(a^(oyo (0, and each one 

is a Rademacher random variable (regardless of the choice of y/^Hi)), which implies ( Gur -L 
(Gt/fll^Mf^M^v) holds unconditionally. Thus in the bootstrapping procedure, we have independence 
among entries within the matrix. 

Let us move to verify the sub-Gaussianity of M matrix. Note that for the index i, j that is not a clique 
for any of the matrices, Mjj is sub-Gaussian, due to Hoeffding’s inequality 

P (|Mij - tMij | > u) < 2exp(-zz 2 /2). (43) 

For the index i, j being a clique in at least one of the matrices, we claim the number of matrices has (z, j) 
being clique is 0*(1). Due to Bernstein’s inequality, we have max; |{0 < s<l: i£ Cjv.kII ^ + § log « 

with probability at least l-zz -1 . This further implies there are at least Z 2 -(^ + | log/r ) 2 many independent 
Rademacher random variables in each z, j position, thus 


P (| Mij - EMij\ >u)< 2exp (-(1 - C- (zczT 1 + l~ l logzz) 2 ) zz 2 /2). (44) 

Up to now we have proved that when z, j is a signal node for M, then O* (1) l~ l > EM,y > l~ l . Thus we 
can take sub-Gaussian parameter to be any o < 1 because kzz -1 , l~ l logzz are both o(l). The constructed 
M{n, k,A/cr) matrix satisfies the submatrix model with XIa = l~ l and sub-Gaussian parameter (7 = 1- 
o(l). 
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Let us estimate the corresponding k in the submatrix model. We need to bound the order of the 
cardinality of R/, denoted as \Ri\. The total number of positions with signal (at least one clique node 
inside) is 


E|i?/| = E|{1 < i < n:i e f?/}| = n 
Thus we have the two sided bound 


1 — (1 -Kin) 


kI 


kI 


1-— <E|2?/| <Kl 
v 2 nl 


which is of the order k:=Kl. Let us provide a high probability bound on \R[\. By Bernstein’s inequality 


P{\\Rl\ — E|i?/|| > u) <2exp 

Thus if we take u = i/4jcZlog«, as long as log n = o(kZ), 

p[\\R l \-E\R l \\>^/4Kl\ogn 


u 2 l 2 
„ kI + uI 3 


(45) 


< 2 n 


(46) 


So with probability at least 1-2 n 1 , the number of positions that contain signal nodes is bounded as 


kI 


kI 

1 - 

v n 


1 4logn 


kI 


< |f?/l <kI 


1 + 


/4logn 


kI 


\Ri\~kI. 


(47) 


Equation (47) implies that with high probability 


kZ( 1 — o(l)) < \R,] <kZ(1 + o(D), 
kZ(1 -0 (D) <\Ci\<Kia + 0(D). 


The above means, in the submatrix parametrization, k m ~ k n ~kI~ n a , XIa ~ Z _1 = n~P, which implies 
k - n a ~P. 

Suppose there exists a polynomial time algorithm that pushes below the computational bound¬ 
ary. In other words, 



(l-2a)/2 . 


B> a — 
2 


(48) 


with the last inequality having a slack e > 0. More precisely, returns two estimated index sets R n 
and C n corresponding to the location of the submatrix (and correct with probability going to 1) under 
the regime /3 = a - 1Z2 + e. Suppose under some conditions, this algorithm s^m can be modified to a 
randomized polynomial time algorithm that correctly identifies the hidden clique nodes with high 
probability. It means in the corresponding hidden clique graph ( S{2N,2k), also pushes below the 
computational boundary of hidden clique by the amount e: 

k{N)=2k~ {2ri) a ~P ~ n ll2 ~ e ^ n 112 - N^. (49) 


In summary, the quantitative computational lower bound implies that if the computational boundary 
for submatrix localization is pushed below by an amount e in the power, the hidden clique boundary is 
correspondingly improved by e. 

Now let us show that any algorithm s4m that localizes the submatrix introduces a randomized al¬ 
gorithm that finds the hidden clique nodes with probability tending to 1. The algorithm relies on the 
following simple lemma. 
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Lemma 13. For the hidden clique model (N, k), suppose an algorithm provides a candidate set S of size 
k that contains the true clique subset exactly. If 

k > C\JklogN 

then by looking at the adjacency matrix restricted to S we can recover the clique subset exactly with high 
probability. 

The proof of Lemma 13 is immediate. If i is a clique node, then min, Jjjec Qy - k-C/2- \J klogN. If 
i is not a clique node, then max, Ljec G,-y < C/2 • yj klogN. The proof is completed. 

Algorithm s4m provides candidate sets Rj, Ci of size k, inside which k are correct clique nodes, and 
thus exact recovery can be completed through Lemma 13 since k £3 (fclogAi ) 1/2 (since k - n l!1 ~ £ £3 k llz - 
n al2 when e is small). The algorithm s£m induces another randomized polynomial time algorithm sdcg 
that solves the hidden clique problem ^(2W,2 k) with k 33 N ] 12 . The algorithm sdcg returns the support 
Cjv.k that coincides with the true support Cn,k with probability going to 1 (a contradiction to the hidden 
clique hypothesis HQ). We conclude that, under the hypothesis, there is no polynomial time algorithm 
.b/m that can push below the computational boundary A 33 

□ 


Proof of Lemma 4. The proof of this lemma uses the well-known Fano’s information inequality, namely 
Lemma 12. We have that X = M+ Z, where M e U mxn is the mean matrix. Under the Gaussian noise with 
parameter a, the probability model is 

P(X|M)oc exp(-<X-M,X-M}/2er 2 ) (50) 

where M = A • UV T e 0. The parameter space 0 is composed of all M = A • UV T where U are sampled 
uniformly on the collection of vectors with k m ones and other coordinates being zero, and similarly V 
are sampled uniformly with k n ones and the rest zero. The cardinality of the parameter space is 


Card(0)= 


/ m 

{ H 

v m , 



corresponding to that many probability measures on the same probability space. Put a uniform prior on 
this parameter space and invoke Fano’s lemma 12. To obtain the lower bound, we need to upper bound 
the Kullback-Leibler divergence ^klO^mII^’) for any Me©, where 


S? - EM'~unif(0)^V- 
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For any Me©, 


^KL^mII^ 8 ) = Esa M log—r- 

o/> 


= E^ m log 

^ ^X~3» u < 

= ^x~e? M i 

1 


&M 


[ (X - M,X - M) 


2 a 2 


1 v {X- M',X — M') 
+ Card( 0 ) J^ 0 


1 


„ (M-M',X-M) 

E -- + 


y 

Card(0) 


Card( 0 ) ]^; £ 0 

{M - M', M - M') 


2 a 2 


2 a 2 I 

1 „ {M— M',M — M') 

Card( 0 ) 


2a 2 


2 a 2 


(M,M) + Car j (Q) <M') ~ 2 Car d(0i 


2 a 2 


A?k m k n k m k n ^ 


a* 


ran 


Thus as long as 


A 2 k m k n k m k n 

(1 -) < a log 


cr / 


mn 


( m 

{ 

\k m 

l kn ). 


we have 


r E ^klI^mII^*) ^ alog(Card(©)). 

Card( 0 ) 


(51) 


(52) 


Invoke the simple bound on binomial coefficients [j) k < (?) < (nr)- If we choose 

A< Ca-cri 


/fc m log£ + fc n log^ 


kin ky 


(53) 


then the condition (28) holds. Any submatrix localization algorithm translates into a multiple testing 
procedure that picks a parameter M' e 0. By Fano’s information inequality 12, the minimax error, which 
is also the localization error, is at least 1 - a - kmlog(m/k ^l\ ilog(n/u • Q 

Proof of Lemma 5. Recall the definition 4 of a sub-Gaussian random variable. Taking Z,y, i e I, j e J, we 
have the following concentration from the Chernoff’s bound for Z/e/j'e/ Ztj 


Ee A Liel.jejZij < exp (|/||/| .(7 2 A 2 /2 c) 


and 


E Zij\>y/\I\\T\(Tt 

iel.jej 


< 2exp(-c- 1 2 12). 


(54) 


(55) 
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There are in total 


km 11 k 


such submatrices, so by a union bound, we have 


m ™ I I Zijl^y/lWlat < -2exp(-c- 1 2 12) 


< 2exp(fc m log(ern/fc m ) + k n log[enlk n ) - ct 2 12). 


If we take t = then with probability at least 


km. /1 k- 


we have 


Thus if 


max | Y, z ij I - ~F v/l !\ I/I °\! log( emlk m ) + k„ log {enlk„ 


\I\=Rm,\J\=C n 


i£l,j£j 


k m k„ A>2 max | Y Z t A (57) 

\i\=km,\J\=k„ i€ f-j €j J 

then the maximum submatrix is unique and is the true one. Recollecting terms, we reach 

4 /log(em/ k in ) log(en//c„) 

- k - + - k -' (58) 

v c v 

To make the proof fully rigorous, we need the following monotonicity trick. Consider the submatrix 
of size k m k n with a rows to be in the correct set R m and b columns to be in the correct set C„ , where 
a < k m and b <k n . The cardinality of the set of such matrices is 


m — a n — a 


km & I 1 kn b 


Using the same calculation as before we want 


0 / {k m - a ) log(e(m - a )/ {k m -a)) + (fc„ - b ) log(e(w - b ) / (/c„ - b)) 

\J~C V kfyi kyi (lb 


By simple algebra, 


kyji d 1 

- < - 

k m k n ab kn 

k n -b 1 

- < — 

k m k n k m 

log(e(m - a) / {k m - a)) < log {em) 
log {e{n-b)l{k n -b)) <log {eri). 
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Hence, if equation (58) is satisfied, (59) is satisfied up to a universal constant for all a < k m and b < k n . 
Thus we have proved that if 


A> C-o\ 


/log {em/k m ) log {enlk n 


y kn k m 

with a suitable constant C, the statistical search algorithm picks out the correct submatrix. The sum of 
the probabilities of the bad events is bounded by 


m- a 

n-b 

k k km ~ Cl 

k n - b) 


-l 


km kn 


(. m-k m )(n-k n ) 


E 

0<a<k m ,0<b<k n 

For the multiple non-overlapping submatrices case, as long as 

m~Y^ k ( s m) - m n - E k^ — n 

5=1 5=1 

then sequential application of Algorithm 4 will find the r-submatrices. 


□ 


Proof of Lemma 3 for Multiple Non-overlapping Submatrices Case. We are going to provide theoretical jus¬ 
tification to the extension of the submatrix localization algorithm to multiple non-overlapping subma¬ 
trices case as in Algorithm 3. Write out the matrix form of the submatrix model, with the SVD version of 
the signal matrix M 

X = M + Z = UAV T + Z. 

Due to the non-overlapping property, we have 1 < 5 ^ t < r, 1 l^ f = 0, so as to C n . The singular values of 

UAV t are k s \j' k { f r> k[ n \ 1 < s< r, and all the other singular values are 0. 

Let us apply the Davis-Kahan-Wedin bound to X = UAV T +Z. Denote the top r left and right singular 
vector of X as U and V. Using Weyl’s interlacing inequality, we have 

|(7 5 (X)-a s (M)|<||Z|| 2 


and 


(7, (X) > A, V k ( s m) k ( f ] - IIZII2, 1 < 5 < r; 
<7r(X) < || Z|| 2 , t> r. 


Thus applying Lemma 7, we have 

max{| sinZ(f/, f7)|,| sinZ(V, V0|} < 


Co{\fm+ \fn) 


min A s y kf n kg - Ca[y/m + \fn) 


fbri) iXn) 

l<s<r 

(j{\fm + yfn) 


min k^/ki^k^ 

1 <s<r 


According to the dehnition of the canonical angles, we have 

max{|| UU T - UU t \\ 2 , II VV T - VV T \\ 2 } < C-■ 


crl\fm + y/ri) 


min k<\/k[ m) k[ n) 

l<s<r 
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Now let us assume we have two observation of X. We use the first observation X to solve for the 
singular vectors U, V, we use the second observation X to project the to the singular vectors U, V. Recall 
X = M + Z = U AV T + Z. Fori <j<n 


W&CjX-j ~ M.j \\ e2 < || & 0 {X.j - M.j ) ||, 2 + || {9> 0 - I)M.j ||, 2 


(64) 


For the first term of (31) because X.j - M.j = Z.j e IR m is an i.i.d. isotropic sub-Gaussian vector, we 
have through Lemma 9, for t = (1 + 1 lc) log n, Z.j e K m , 1 < j < n and r > 0 


||^ & (Xy-M 7 -)||, 2 >av'r 




i + 2vTTT7 c - 


/log n , logn 

——+ 2(l + l/c)-—— 


< n 


-c-l 


(65) 


Thus invoke the union bound for all 1 < j < n 


max \\£Pjj(X. j - M.j)\\e 2 < cr\fr + Vl + 1 lc-a\ /log n + (1 + 1/c) 

1 <i<n * 


log 71 

7 ^ 


< a\fr + C • er \/log n 


( 66 ) 

(67) 


with probability at least 1 - n~ c . 

For the second term M.j = X.j - Z.j of (31). It can be estimated through perturbation Sin Theta 
Theorem 7. Basically it is 


|| {9> 0 -1)M.j || h = || - &>u)M.j \\ e2 < || UU T - UU T || 2 • || My || , 2 


<C- 


a\/m+ n 


min A s v/ k ( s k s 


1 <s<r 


max A ,\/ k[ m \ 


Combining all the above, we have with probability at least 1 — n c - m c , for all 1 < j < n 


W 0 x.j-M.j\\ e2 <c- 

Similarly we have for all 1 <i<m, 


\\&vXl-Ml\\e 2 <C- 


o \[r + a \/log n + a\fm\rn ■ 


max A s \/ k[ m) 

t<s<r v * 


min Asx/k^k^ 

1 <s<r v 


o\fr + o\ \ogm + o\fm\rn- 


max A,\/ k['^ 

1 <s<r 


min A s \/ k[ mi k { f ] 

l<s<r 


Clearly we know for any 1 < s < r and i e R s and i' e [m]\R S 


( 68 ) 

(69) 


(70) 


(71) 


I M [. - M,v. ||, 2 > min A S \J k [ s n) 


1 <s<r 


and for any 1 < s < r and j e C s and j' e [n\\C S 


\M.i- M.;i\\f,> min A s \/ k[ m) . 

J J l<s<r v 
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Thus if k[ m ^ - k m , k^ ~ k n and A 5 - A for all 1 < 5 < r 


^n/^ot ^ 6C- cr\fr + o J\ogn + cr 


' m v n 


h\[k~n > 6C- 


cr\fr + 0 \ \ogm + cr\ 


1 mv n 
km 


(72) 

(73) 


We have learned a metric d (of intrinsic dimension r) such that under this metric, data forms into clusters 
in the sense that 

2 max | dj - df \ < min | d t - di' \. 
i,i'cR s ie.R m ,i'e[m]\R s 

Thus it satisfies the geometric separation property. 

Thus in summary if 


A> C-cr 


km A kn 


■ + ■ 


/logn /log m mv n 


kn 


kn 


+ ■ 


ktnkn 


the spectral algorithm succeeds with probability at least 

1 - m~ c - n~ c - 2exp (-c(m + n)). 


Due to the fact that 


1 mv n 


km A k n V k m k n 

because rk m ^ m, rk n ^ n in most cases, the first term does not have an effect in. most cases. 


□ 


Proof of Theorem 3 is a direct result of Lemma 1 and Theorem 2. Proof of Theorem 4 is obvious based 
on Lemma 2 and the hidden clique hypothesis HC|. Proof of Theorem 5 combines the result of Lemma 5 
and Lemma 4. 


References 

Arias-Castro, E., Candes, E. J., Durand, A., et al. (2011). Detection of an anomalous cluster in a network. 
The Annals of Statistics, 39(l):278-304. 

Balakrishnan, S., Kolar, M., Rinaldo, A., Singh, A., and Wasserman, L. (2011). Statistical and compu¬ 
tational tradeoffs in biclustering. In NIPS 2011 Workshop on Computational Trade-offs in Statistical 
Learning. 

Bennett, G. (1962). Probability inequalities for the sum of independent random variables. Journal of the 
American Statistical Association, 57(297):33-45. 

Berthet, Q. and Rigollet, P. (2013a). Computational lower bounds for sparse pea. arXiv preprint 
arXiv:1304.0828. 

Berthet, Q. and Rigollet, P. (2013b). Optimal detection of sparse principal components in high dimension. 
The Annals of Statistics, 41(4):1780-1815. 


29 
























Birnbaum, A., Johnstone, I. M., Nadler, B., and Paul, D. (2013). Minimax bounds for sparse pea with noisy 
high-dimensional data. Annals of statistics, 41 (3): 1055. 

Butucea, C. and Ingster, Y. I. (2013). Detection of a sparse submatrix of a high-dimensional noisy matrix. 
Bernoulli, 19(5B):2652-2688. 

Butucea, C., Ingster, Y. I., and Suslina, I. (2013). Sharp variable selection of a sparse submatrix in a high¬ 
dimensional noisy matrix. arXiv preprint arXiv:1303.5647. 

Cai, T. T., Liang, T., and Rakhlin, A. (2014). Geometrizing local rates of convergence for linear inverse 
problems. arXiv preprint arXiv:1404.4408. 

Cai, T. T., Ma, Z., and Wu, Y. (2013). Sparse pea: Optimal rates and adaptive estimation. The Annals of 
Statistics, 41(6):3074-3110. 

Cai, T. T., Ma, Z., and Wu, Y. (2015). Optimal estimation and rank detection for sparse spiked covariance 
matrices. Probability Theoiy and Related Fields, page to appear. 

Candes, E. J. and Plan, Y. (2011). Tight oracle inequalities for low-rank matrix recovery from a minimal 
number of noisy random measurements. Information Theory, IEEE Transactions on, 57(4):2342-2359. 

Chandrasekaran, V. and Jordan, M. I. (2013). Computational and statistical tradeoffs via convex relax¬ 
ation. Proceedings of the National Academy of Sciences, 110(13) :E1181-El 190. 

Chandrasekaran, V, Recht, B., Parrilo, P. A., and Willsky, A. S. (2012). The convex geometry of linear 
inverse problems. Foundations of Computational Mathematics, 12(6):805-849. 

Chen, Y. and Xu, J. (2014). Statistical-computational tradeoffs in planted problems and submatrix local¬ 
ization with a growing number of clusters and submatrices. arXiv preprint arXiv:1402.1267. 

Deshpande, Y. and Montanari, A. (2013). Finding hidden cliques of size\ sqrt {N/e} in nearly linear time. 
arXiv preprint arXiv:1304.7047. 

Donoho, D. and Jin, J. (2004). Higher criticism for detecting sparse heterogeneous mixtures. Annals of 
Statistics, pages 962-994. 

Donoho, D. L. and Johnstone, I. M. (1998). Minimax estimation via wavelet shrinkage. The Annals of 
Statistics, 26(3):879-921. 

Feldman, V, Grigorescu, E., Reyzin, L., Vempala, S., and Xiao, Y. (2013). Statistical algorithms and a 
lower bound for detecting planted cliques. In Proceedings of the forty-fifth annual ACM symposium on 
Theory of computing, pages 655-664. ACM. 

Gross, D. (2011). Recovering low-rank matrices from few coefficients in any basis. Information Theory, 
IEEE Transactions on, 57(3): 1548-1566. 

Hoeffding, W. (1963). Probability inequalities for sums of bounded random variables. Journal of the 
American statistical association, 58(301):13-30. 

Hsu, D., Kakade, S. M., and Zhang, T. (2012). A tail inequality for quadratic forms of subgaussian random 
vectors. Electron. Commun. Probab, 17(6). 


30 



Kolar, M., Balakrishnan, S., Rinaldo, A., and Singh, A. (2011). Minimax localization of structural informa¬ 
tion in large noisy matrices. In Advances in Neural Information Processing Systems, pages 909-917. 

Latala, R. (2005). Some estimates of norms of random matrices. Proceedings of the American Mathemat¬ 
ical Society, 133(5):1273-1282. 

Lee, D. D. and Seung, H. S. (2001). Algorithms for non-negative matrix factorization. In Advances in 
neural information processing systems, pages 556-562. 

Ma, Z. and Wu, Y. (2013). Computational barriers in minimax submatrix detection. arXiv preprint 
arXiv:1309.5914. 

McSherry, F. (2001). Spectral partitioning of random graphs. In Foundations of Computer Science, 2001. 
Proceedings. 42nd IEEE Symposium on, pages 529-537. IEEE. 

Montanari, A. and Richard, E. (2014). Non-negative principal component analysis: Message passing 
algorithms and sharp asymptotics. arXiv preprint arXiv:1406.4775. 

Ng, A. Y., Jordan, M. I., Weiss, Y., et al. (2002). On spectral clustering: Analysis and an algorithm. Advances 
in neural information processing systems, 2:849-856. 

Shabalin, A. A., Weigman, V J., Perou, C. M., and Nobel, A. B. (2009). Finding large average submatrices 
in high dimensional data. The Annals of Applied Statistics, pages 985-1012. 

Stewart, G. W. and Sun, J.-g. (1990). Matrix perturbation theory. Academic press. 

Tsybakov, A. B. (2009). Introduction to nonparametric estimation, volume 11. Springer Series in Statistics. 

Tsybakov, A. B. (2014). Aggregation and minimax optimality in high-dimensional estimation. 

Vershynin, R. (2010). Introduction to the non-asymptotic analysis of random matrices. arXiv preprint 
arXiv:1011.3027. 

Vu, V (2008). Random discrete matrices. In Horizons of combinatorics, pages 257-280. Springer. 

Vu, V (2014). A simple svd algorithm for finding hidden partitions. arXiv preprint arXiv:1404.3918. 

Vu, V. Q. and Lei, J. (2012). Minimax rates of estimation for sparse pea in high dimensions. arXiv preprint 
arXiv:1202.0786. 

Wainwright, M. J. (2014). Structured regularizers for high-dimensional problems: Statistical and compu¬ 
tational issues. Annual Review of Statistics and Its Application, 1:233-253. 

Wang, T., Berthet, Q., and Samworth, R. J. (2014). Statistical and computational trade-offs in estimation 
of sparse principal components. The Annals of Statistics, to appear. 

Zass, R. and Shashua, A. (2006). Nonnegative sparse pea. In Advances in Neural Information Processing 
Systems, pages 1561-1568. 

Zhang, Y., Wainwright, M. J., and Jordan, M. I. (2014). Lower bounds on the performance of polynomial¬ 
time algorithms for sparse linear regression. arXiv preprint arXiv:1402.1918. 


31 



A Convex Relaxation Algorithm 


In this section we will investigate a convex relaxation approach to the problem. The same algorithm 
has also been investigated in a parallel work of Chen and Xu (2014). Our analysis is slightly different, 
with the explicit construction of the dual certificate using the idea in Gross (2011). For the purposes of 
comparing to the spectral approach, we include the convex relaxation analysis in this section. Let us 
write the optimization problem 


min ||X- Auv T \\ % 

weR^i/eR" r 

s.t. \\u\\e Q = k m ,\\v\\e a = k n 
u £ {0 ,\} m , v e {0,1}". 

This problem is non-convex: the feasibility set is non-convex, and so is the optimization function (al¬ 
though it is bi-convex). However, we can relax the problem and transform it into a convex optimization 
problem. Of course, we need to ensure that the solution to the relaxed problem is the exact solution (with 
high probability) under appropriate conditions. 

The matrix version of the submatrix problem suggests that the signal matrix is of the structure “low 
rank and sparsity on the singular vectors.” We recall from the low rank matrix recovery literature, see e.g. 
Candes and Plan (2011) and Cai et al. (2014), that we can utilize the low rank structure and solve relaxed 
versions as follows. 


Relaxation 1 Consider the constraint minimization relaxation, 


min 

MeR" 1 *" 

S.t. 


IIMIU 

\\X-XM\\ 2 <C-o{s/m + y/TD. 


Unfortunately, Relaxation 1 is only good in terms of estimation of the whole matrix. The stronger ob¬ 
jective of localization requires simultaneous exploitation of sparsity and low rank-ness, as in Relaxation 

2 . 


Relaxation 2 Let us expand the objective of the original non-convex optimization problem, drop the 
quadratic term to make the procedure adaptive in terms of A, and convexify the feasibility set at the 
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same time. 


Algorithm 5: Convex Relaxation Algorithm 


Input: X £ R' nxn the data matrix. Size of the submatrix k m x k„. 

Output: A subset of the row indexes R m and a subset of column indices C n as the localization sets 


of the submatrix. 

1. Solve the following convex optimization problem 


Mq = argmin - (X, M) 


MeR mx " 


s.t. 0<M<l m ll 


l|M||* < 1 

(M, l m l n ) = k m k n 


2. Perform SVD on Mq, denote the set of the non-zero entries on the top left singular vector to be 


R m and the set of the non-zero entries on the top right singular vector to be C n . 


The time complexity to solve this convex optimization problem is at least 0((m + n) 3 ) implemented 
with alternating direction methods of multipliers (ADMM). The disadvantage is that the theoretical guar¬ 
antee only holds for the exact solution Xy, however, in reality we can only approximately find X (l through 
ADMM or some other optimization methods. We also remark that this algorithm requires the prior 
knowledge of the submatrix size k m , k n , which means it is not fully adaptive. 

Lemma 14 (Guarantee for Relaxation Algorithm). Consider the submatrix model (2) and the Algorithm 5. 
There exists a universal C > 0 such that when 



the convex relaxation succeeds (in the sense that R m = R m , C n = C n ) with probability at least 1 - 2m c 


2 n c -2 [mn] e - 2exp (-c(m + n)). 

Proof of Lemma 14. Let us construct the dual certificate to secure that the true solution is the unique so¬ 
lution. Ifwe can construct a pair of a primal certificate M* = \J k m k n UV T and dual certificate A*, 0*, p*, v* 
that satisfy 



(74) 

(75) 


(76) 


Here W £ SP U L nV i, || W H 2 < 1, and UV T + W denotes the sub differential of || • II * evaluated at M* 


d\ m=m* IIMII* = UV t + W. 


Equation (75) is equivalent to 


A *j > 0 , e*j = 0, it Rm U jtCn 
&■ j > 0, A- j = 0, i £ R m n j £ C n . 
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(77) 

(78) 













We claim that any solution M to Relaxation 2 must satisfy M = M *. If not, we write M = M* + H and find 
that 


0 < M* + H< l m l% (79) 

|| M* + H || * < || M* || * => (H, UV T +W)<0 (80) 

{H,l m \ T n ) = 0. (81) 

All of the above equations are due to the primal feasibility, and the second inequality also uses the con¬ 
vexity of || -1| *. Note that (79) can be written in a more explicit form 

0 < Hij <1, it R m u jtC n (82) 

-1 < Hij <0, i e R m n je C n . (83) 

Due to the optimality of M in terms of the objective function 

(X, M* + H) > (X, M* > 


which means 


0 <(X,H) (84) 

= <—A* + 0* + /d* [UV T + W) + v* l m ll, H ) (85) 

<<—A* + Q* ,H). ( 86 ) 


We can see that if Hjj > 0, then we must have M*. = 0, which through complimentary slackness implies 
0f. = 0. This in turn means - A/ ,■ < 0. When Hj / < 0, we must have M*. = 1, which again means A,, = 0, 

IJ J J l J J 

0*, > 0, a contradiction. Thus Hj , = 0 for all i, j. 

IJ J 

The properties we impose on the dual certificates are motivated from the Karush-Kuhn-Tucker (KKT) 
conditions. Introduce the dual variables A,0 e IR" !X ”, p e K + , v e R for the four feasibility conditions. 
Then the Lagrangian is 

^(M,A,0,p,v) = -(X,M)-(A,M> + (0,M-l m lJ> + p(||M|U-l) + v((M,l m lJ>-fc m fc„). 

The associated KKT conditions are 


X — — A + 0 + ii{UmVm + Wm) + (87) 

0<M<l m ll (88) 

I|M|U<1 (89) 

(ikf, \m^ n ) = k m kn (90) 

A>0,©>0,p>0 (91) 

A ijMij = 0,&ij (1 - Mij ) = 0. (92) 


Now let us see how to construct the dual certificate - A* + 0*, p*, v* to satisfy the conditions (74) - 
(76). Expand (87) as 


-A* + 0* = \l Rm \ T Cn + &UU V (Z) + &uT nV T (Z) - p* 


\J k m k n 


l« m l C n +W 




(93) 
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where 


& a r nV T{Z) = Z-& UuT {Z). 


(94) 


Choose 


W=—& u r nV T(Z). 


(95) 


Thus if we choose p* > C- o{\fm + \fn), it holds that 

\\W\\2 = ^\\&>u T nV T (Z)\\2<^\\Z\\ 2 <l 
ji* p* 

with probability at least 1-2 exp (-c(m + n)). Thus with the choice of W, Equation (93) becomes 


A*+0*=Al« m l' n + 5W {Z)-n* 


1 


\J k m k, 




Hence, we need to have 


A-p* 


_ -v* -max\[@ > UuV {Z)] i j\ >0, i e R m n jeC„ 

V k m kn ij 

max\[3? UuV {Z)]ij\-v* <0, i e R c m u jeC c n . 
ij 


Now let us write out the explicit form of the projection S^uuv(Z) 

1 


&Uuv(Z) = -Li Rm il m Z+ Z\ Rn \ T Rn - —— 
Let us see the concentration property of [£Puuv(Z)]ij: 


l* m l Lzin.lL 


(96) 

(97) 

(98) 

(99) 


mjuviZftij = 

\[£? > u\jv(Z)]ij \ < 

For all the 1 < j < n 


An 


E Zkj 

\k€.R m j 


1 icRm + ' 


E zu 

\lcCn } 


j£-C n 


E F 


kmkn [keR^leC,, 1 1 


i€.R m >j E-Cn 


f E z ^ 

Km kER m 


+ 


T E zu 

IeC„ 


+ 


1 


- E Zu 

K mKn keR, n ,leC n 


( 100 ) 

( 101 ) 


max 

l<j<n 


f E 

Km kER m 


< \/2(l + l lc)- 


(J\ 


llogn 


( 102 ) 


with probability at least 1 - 2n c . For all the 1 < i < m 


max 

1 <i<m 


F E 

Kn IeC„ 


< v/2(l + l/c)-o-i 


/logm 


(103) 


with probability at least 1 - 2 m c . For all i, j, 

1 


max 

ij 


E Zij 


kmkn kER m ,lEC n 


< V2(l + 1 !c)-o\ 


/log(mn) 
km kn 


(104) 
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with probability at least 1 - 2{mri) c . Now, pick the dual certificate variables in the following way 


p* = C-o{\fm + \fn) (105) 

... „ /log"* , / lo 8 n \ , inr , 

v =CCT lV^r V^rJ H06) 

0* >0, ie R m n jeCn (107) 

A* <0, ie R c m u jeC c n (108) 


where the last two equation follows from (97) and (98). We conclude that the relaxation algorithm suc¬ 
ceeds with probability at least 

1 -2 m~ c - 2 rT c -2 (mn)~ c - 2 exp(-c(m + n)) 


if 


A> C-o 


t 


l m + n 

km ki 7 


■ +1 


logm /log n 


kti 


■ + 


kt i 


We have achieved the the same boundary as the spectral method upper bound. 


□ 


B Algorithmic Reduction for Detection 


Theorem 6 (Computational Lower Bounds for Detection). Consider the submatrix model (2) with pa¬ 
rameter tuple (m = n, k m — k n ~ n a ,X/a = n~k), where / < a < 1, /)>(). Under the hardness assumption 
HQ, if 


A . m + n 


it is not possible to detect the true support of the submatrix with probability going to 1 for any polynomial 
algorithm. 


Proof of Theorem 6. We would like to build a randomized polynomial mapping from the hidden clique 
graph < £{N,k(N)) to a matrix M{m = n, k m — k n - k, XIo) for the submatrix model. Denote this transfor¬ 
mation as 

ST : <0{N,k{N)) — M{n, k = n a ,X/a = n~h- 

There are several stages of the construction. First, we define a graph that is stochastically equiva¬ 
lent to the hidden clique graph C S, but is easier for the analysis. Let us call it c S e . c S e has the property: 
each node independently has the probability k{N)I N to be a clique node. By Bernstein’s inequality, with 
probability at least 1 - 2IV -1 , the number of cliques K e in c S e 


r 

k 1 


I 4log N 


<K <K 


1 + - 


/ 4logAi 


(109) 


as long as k £ log N. 

Consider a double sized hidden clique graph ( S e [2N,2 k{N)) with N = n 1+ P, and k(N) = k = n a , \ < 
a < 1. Denote the clique nodes set as Cn, k - Connect the hidden clique graph to form a symmetric matrix 
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G £ {-1, i} 2JVx2JV , where Gjj = 1 if i,j e Cn,k> otherwise with equal probability to be either -1 or 1. Take 
out the upper-right submatrix of G, Gur where U is the index set 1 < i < N and R is the index set N + 1 < 
3 * 2 N. 

Partition the rows of Gur £ R” xn , to form n blocks. The s’s block, 1 < s < n, corresponds to the 
row index set I s = {i : (s- l)nP + 1 < i < srJ’\. Construct the M e U nxn matrix in the following way 

S'- M st = — £ (GuR)ij, 1 <s,t<n. (110) 

nP i£l si j£l, 

There are two cases, if I s n C/v, K # 0 and I t n Cn,k 0 0. namely when s’s block contains at least 1 clique 
node, and so does t’s block, then E M st > n~P. Otherwise E M st = 0. Note that for the index s, t that has no 
clique inside, E M st is sub-Gaussian random variable, due to Hoeffding’s inequality 

P{\M st -EM st \ > u) <2exp(-w 2 /2). (Ill) 


For the s, t with clique nodes inside, we know the maximum number of clique nodes is logn (due to 
Bernstein’s inequality that maxi< 5 < n \ I s n C,\! iK | < -^ +1 log « with probability at least 1 - n ~ 1 ), which means 
there are at least ^ -1 log n many independent Rademacher random variables in each s, t block, thus 

P{\M sr -EM sr \ > u) <2exp^-(l-C-(kn _1_ ^ + n _ ^logn))u 2 /2j. (112) 

Thus we can take sub-Gaussian parameter to be any cr < 1 because kn~ l ~P , n~^\ogn are both o(l). Now 
this constructed M{n, k ) matrix satisfies the submatrix model with A = n~P and sub-Gaussian parameter 
£7=1- o(l). 

Let us see how many elements in 1 < 5 < n are such that I s n Cn,k 0 0- Namely, we want to estimate 
how many clique nodes there exist in the transformed submatrix model. We have the two sided bound 


k 



1 

2n l ~ a , 


< E| {s: I s n C n , k ^0,l<s<n}\<k. 


which is of the order k. Using Bernstein’s bound, we have with high probability 


1- 


1 ) 


2 n l ~ a ) 


1- 


/4logn 




< \{s:I s rC N ,K 7^0,l<s<n}\<k 


1 + ' 


1 4logn 


Thus the submatrix model M(m = n, k m - k n , XIa) satishes k m - k n ~ k. 

Suppose there exists an polynomial time algorithm s^m that pushes below the computational bound¬ 
ary quantitatively by a small e > 0 amount 


r~P x - 


m + n 


,l-2a 


/3 > 2a - 1. 


(113) 


where /3 = 2a - 1 + e. Namely, &4m detects below the boundary then it naturally introduced a 

polynomial time detection algorithm for hidden clique problem < £(N= n l+ P, k = n a ) , which violates the 
HQ because 

k{N) ~ N = IV^ NK 

□ 
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